#%%
"""
Created on Jan 22 2019
Monte Carlo paths for the Jacobi correlation model
@author: Lech A. Grzelak
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d 


def mainCalculation():
    t1=[0,0.821900,1.64380,2.46580,3.28770,4.10960,4.93150,5.75340,6.57530,\
        7.39730,8.21920,9.04110,9.86300,10.68490,11.50680,12.32880,13.15070\
        ,13.97260,14.79450,15.61640,16.43840,17.26030,18.08220,18.90410,19.7260,\
        20.54790,21.36990,22.19180,23.01370,23.83560,24.65750,25.47950,26.30140,\
        27.12330,27.94520,28.76710,29.5890,30.4110,31.23290,32.05480,32.87670,\
        33.69860,34.52050,35.34250,36.16440,36.98630,37.80820,38.63010,39.45210,\
        40.2740,41.09590,41.91780,42.73970,43.56160,44.38360,45.20550,46.02740]
    t2=[0,0.821900,1.64380,2.46580,3.28770,4.10960,4.93150,5.75340,6.57530,7.39730,\
        8.21920,9.04110,9.86300,10.68490,11.50680,12.32880,13.15070,13.97260,\
        14.79450,15.61640,16.43840,17.26030,18.08220,18.90410,19.7260,20.54790,\
        21.36990,22.19180,23.01370,23.83560,24.65750,25.47950,26.30140,27.12330,\
        27.94520,28.76710,29.5890,30.4110,31.23290,32.05480,32.87670,33.69860,\
        34.52050,35.34250,36.16440,36.98630,37.80820,38.63010,39.45210,40.2740,\
        41.09590,41.91780,42.73970,43.56160,44.38360,45.20550,46.02740]
    df1= [1,0.9912650,0.9805340,0.9641250,0.9418240,0.9192670,0.895380,0.8692230,\
          0.8428550,0.8164780,0.7899110,0.7632080,0.7367020,0.7107840,0.684920,\
          0.6594810,0.6349270,0.6105070,0.5862780,0.5650440,0.545270,0.5258950,\
          0.5069240,0.4883650,0.4702220,0.4546480,0.4408020,0.4273960,0.4144180,\
          0.4018510,0.3896830,0.3789860,0.3695590,0.3605050,0.3518070,0.343450,\
          0.3354190,0.3282840,0.3221350,0.3162880,0.310730,0.3054490,0.3004350,\
          0.2956780,0.2911670,0.2868940,0.282850,0.2790270,0.2754190,0.2711330,\
          0.2648390,0.2587240,0.2527820,0.2470070,0.2413950,0.2359390,0.2306360]
    df2=[1,0.9983970,0.9968090,0.994460,0.9907650,0.985370,0.977910,0.9684380,\
         0.9567750,0.9430350,0.927610,0.9110590,0.8937890,0.8760590,0.858030,\
         0.8402370,0.8225770,0.8051410,0.7881340,0.7716530,0.7555510,0.7398430,\
         0.724550,0.7096920,0.6952860,0.6813370,0.6678080,0.6546720,0.6419040,\
         0.6294840,0.617390,0.6056190,0.5942220,0.583150,0.5723440,0.5617510,\
         0.5513190,0.5410060,0.530820,0.5207770,0.5108910,0.5011740,0.4916370,\
         0.4822890,0.4731410,0.4641990,0.4554720,0.4469650,0.4386860,0.4306370,\
         0.4228050,0.4151740,0.407730,0.4004620,0.3933560,0.3864040,0.3795930]

    plt.figure(1)
    plt.plot(t1,df1);
    plt.plot(t2,df2,'--r');
    plt.grid()
    plt.xlabel('T [years]')
    plt.ylabel('P(0,T)')
    plt.title('Zero-Coupon Bond, P(0,T)')
    plt.legend(['P(0,T), mid 2009','P(0,T), mid 2014'])
    plt.axis([0.0,np.max(t1),np.min(df1),np.max(df1)])

    # Shock size

    dt     = 0.01

    # Define ZCBs as the interpolated function of market data

    P_0_t1 = interp1d(t1,df1,fill_value='extrapolate') 
    P_0_t2 = interp1d(t2,df2,fill_value='extrapolate')

    print(P_0_t1(4.0))
    print(P_0_t2(4.0))

    # Instantaneous forward rate f(0,t)

    f_t1   = lambda t:-(np.log(P_0_t1(t+dt))-np.log(P_0_t1(t-dt)))/(2.0*dt)
    f_t2   = lambda t:-(np.log(P_0_t2(t+dt))-np.log(P_0_t2(t-dt)))/(2.0*dt)

    plt.figure(2)
    plt.plot(t1,f_t1(np.array(t1)))
    plt.plot(t2,f_t2(np.array(t2)),'--r')
    plt.grid()
    plt.xlabel('T [years]')
    plt.ylabel('f(0,T)')
    plt.title(['Instantaneous forward rate, f(0,T)'])
    plt.legend('f(0,T), mid 2009','f(0,T), mid 2014')
     
mainCalculation()
